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ABSTRACT 

A Lunar Gradiometer mission Involves the mounting of 
a gradiometer on a satellite which is in a low, polar, 
and circular lunar orbit. This paper presents the re- 
sults of a numerical simulation of the mission. It is 
shown that if the satellite is in a 50 km orbit, 1° and 2° 
gravity anomalies may be estimated with accuracies of 
12 mgal and 1 mgal respectively. At a 100 km altitude, 

2° gravity anomalies can be estimated with an accuracy 
of 12 mgal. These results assume a rotating type gradiom- 
eter with a . IE accuracy. The results can be readily 
scaled to reflect another accuracy level. 
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SIMULATION OF A LUNAR GRADIOMETER MISSION 


INTRODUCTION 

For the purpose of estimating the lunar gravity field, the output of a padiometer 
has advantages over satellite perturbation data. The gradiometer, since it can 
measure second derivatives of the lunarpotential field, will be sensitive to the 
density variations of the Moon’s outer crust which are the source of the signifi- 
cant high frequency components of the lunarpotential field. These high frequency 
effects tend to be smoothed out in satellite perturbation data. 

Regardless of what data type one employs there are difficulties in using standard 
parameter estimation techniques to determine a global and detailed lunar gravi y 
field. The essence of these difficulties is that a large number of parameters 
must be estimated. For instance, if features of the lunar potential as small as 
2“ are to be recovered and if the standard spherical harmonic expansion of the 
lunarpotential field is used as a parameterization, a full set of coefficients to 
degree and order 90 is required. This implies the estimation of over 8 000 param- 
eters. In general it is not possible to simultaneously estimate such large param- 
eter sets. In practice, in order to use standard estimation techniques to recover 
lunarpotential fine structure from gradiometer data, it will be necessary to ad- 
just small subsets of parameters while constraining other parameters to a priori 
values. But unless the parameterization exhibits a certain orthogonality property 
in the data set the effect is that the uncertainties in the unadjusted terms will 
badly corrupt the estimates of the adjusted terms. If spherical harmonic coef- 
ficients are utilized the estimation problem can be decomposed into estimation 
problems of smaller dimensionality only if the gradiometer data set is globaUy 
distributed. This implies a heavy computational load. A parameterization of 
the lunarpotential field is desired in which the individual parameters exhibit a 
localized observability pattern in gradiometer data. With such a parameteriza- 
tion it would be possible to accurately estimate small subsets of parameters by 
processing localized blocks of gradiometer data. The mean free air gravity 
anomaly - Stokes’ function parameterization of a gravity field [1] has been shown 
to possess a measure of orthogonality in local blocks of altimeter data [2]. In 
[3] the observability of gravity anomalies in gradiometer data was shown but the 
orthogonality property was assumed rather than demonstrated. To be realistic 
one must recognize that the orthogonality is not perfect and that only the esti- 
mates of gravity anomalies which are a sufficient distance from unadjusted grav- 
ity anomalies will be of value. This implies that gravity anomalies should be 
estimated in blocks with the estimates of anomalies in a sufficiently small inner 
core accepted as valid and the rest rejected due to aliasing. To obtain an intelli- 
gent estimation strategy for estimating lunar gravity anomalies from gradiometer 
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data it is necessary to determine for a given data block size the relationship be- 
tween the accuracy with which a given anomaly is estimated and its distance from 
the nearest imadjusted anomaly. 

In this report lunar gradiometer missions at 50 km and 100 km altitude are sim- 
ulated and optimal data block and estimation block sizes for determining 1” a nr? 

2“ lunar gravity anomalies are obtained by means of covariance analysis* A 
rotating type gradiometer as described in [4] is assumed mounted on a lunar 
polar satellite and to have an accuracy of . IE, (IE = 10~^ gal/cm). The results 
can be readily scaled to reflect any other accuracy level. 
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ORTHOGONALITY PROPERTIES OF GRAVITY ANOMALIES IN 
GRADIOMETER DATA 

The rotating gradiometer is an. instrument which senses gravity gradients in a 
given plane. If the instrument is on board a satellite we take that plane to be the 
orbital plane of the satellite. Define a local satellite centered coordinate system 
(Ip Ij, I 3 ) where the Ij unit vector is pointii^ northward, the Ij unit vector is 
pointing eastward, and the 1 3 unit vector is perpendicular to the plane spaced by 
I j gnd 1 2 and is directed outward from the moon. If the satellite on which the 
rotating gradiometer is mounted is in a polar orbit, the sensing plane of the in- 
strument always coincides with the plane spanned by the 1^ and Ij unit vectors. 
This assumption considerably simplifies the following mathematical development. 
The lunarpotential field can be represented as 

W = U + T (1) 


where U is a reference lunarpotential generally defined by a low degree and order 
spherical harmonic expansion of a nominal field and T is the so-called anomalous 
potential. The anomalous potential represents the high frequency part of the 
field and we are primarily interested in the portion of the gradiometer output 
which is due to T. Assuming that the satellite is in a polar orbit, the additive 
component of the rotating gradiometer signal amplitude which is due to the anom- 
alous potential T is 


AMP 


/ p^T(r,0A) _I_ d^T(r,^,X)y 3^ T(r,0,X)y V'^ 

\ _ dr^ d<l>^ _ _r 3r90 _ / 


( 2 ) 


where r, <j>, and X have the usual spherical coordinate interpretations and where 
the derivatives are evaluated at the satellite position [3]. 

The anomalous potential T at any point above the surface of the moon can be ex- 
pressed by means of the discrete form of Stokes’ formula as 

' T(r, X) = X') S(r, 0 , X, <(>\, cos ( 0 |) A<.| AX' (3) 


where 5 gj (<{>[, Xj ) is a mean gravity anomaly over a block centered at latitude 
<!>[; and longitude X| and referenced to a equipotential surface defined by the 
nominal field U. The expression cos A(j>[ AX' represents the area of the 
block on which the i‘*^ gravity anomaly is defined and S(r, (j>, X, 0 ?, X|) is Stokes' 
function given in [ 3 ] by 


S(r, 0, X, <j)\ X') = t 


A + 1 - 3 D - t cos |5 + 3?n 


1 + 

_ 


t cos (i//) + D 


(4) 
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where 


i|/ = cos"^ [sin (0) sin (0*) + cos (0) cos (0') cos (V - X)] 


t = — 
r 


D = (1 - 2 t cos(0) + t^)'^ 


where R is the mean radius of the moon. The functional relationship between 
the anomalous part of the gradiometer output and a globally distributed set of 
gravity anomalies is obtained from Equations 2 and 3 as 



' B^S(r,0.X,0;,Xp 

3r2 


1 3^S(r,0,X.0|,Xp '[^ ^ 3^S(r,0,X,0?,X;) P\ 

30^ . 3r30 _ / 


V4 

(5) 


cos (0() A0| AXj 5g(0|,Xj) 


Equation 5 shows that if the discrete form of Stokes' formula is valid, the gradi- 
ometer output is a linear function of a global distribution of gravity anomalies. 
Hence standard parameter estimation techniques should be applicable to the 
problem of determining gravity anomalies from gradiometer data. To see how 
these techniques are applied it is useful to employ matrix notation. Let y be a 
vector of the anomalous components of gradiometer readings. Assume that the 
i^^ component ofj^ is the anomalous gradiometer reading obtained when the satel- 
lite's position is given in spherical coordinates as (rj, 0j, Xj). Let'lf be a vector 
of numerical values of a global set of gravity anomalies. The component of 
g is the numerical value of the gravity anomaly centered at latitude 0j and longi- 
tude Xj . The functional relationship between y and g is expressed as 

y = Ag (6) 


where A is a matrix the number of whose rows is the number of observations 
and the number of whose columns is the number of gravity anomalies. The ele- 
ment in the i^*^ row and colunrn of A is 


A(I,J) = 


~ 3^S(r,0j,?4,0j,x;) 2 8^S(r^,0i,X^,0;,Xj)' 

r? 30? 

cos(0j) A0j AXj 


3r? 


i>T r* 

L +4 _ 

J L^-i 


a2s(ri,0i,Xi,0:Xj)' 

arj 30J 


n'Y 

- / (7) 


Equation 7 provides a linear equation of condition and in a standard minimum 
variance fashion 'g could be estimated from observations of y . In order for Equa- 
tion 7 to be correct (correct, that is, assuming that the approximations inherent 
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in the discrete form of Stokes' formula are valid) the gravity anomalies in the 
array g must cover the moon’s surface. Computational considerations constrain 
us to choose a region such that the number of gravity anomalies to be estimated 
does not exceed a few hundred. In effect gravity anomalies outside of this re- 
gion are assumed to be zero. To see precisely what happens when this assump- 
tion is made postulate the g array of Equation 6 to be defined over the globe, 
and write 



where 

■g^ = gravity anomalies to be adjusted in a standard minimum variance 
filter and 

gj == gravity anomalies assumed to be zero and thus left unadjusted by the 
minimum variance filter. 

Then Equation 6 can be written 


y 




( 8 ) 


where Aj and Aj are respectively the variational matrices of y with respect to 
gi and g" 2 . 


The gradiometer output provides direct observations y of y with statistics 

y = y + E (»/) = 0, E (t" = Q (9) 


The gravity anomalies g'^ for computational reasons are assumed to be zero but 
the actual values of gravity anomalies in the region of the sphere which is ignored 
have a certain distribution about zero. We assume 

E (g^ ) = 0, E (gj gj) == Pj (10) 


When the values of gj are assumed to be zero the minimum variance estimate 
of gj becomes 

t, = (ATq-Ia^)-i A^Q-iy (11) 

Define the covariance matrix of the estimator given by Equation 11 as 

P = E(l^, -gj] [^j -gj]T) (12) 
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From Equations 9, 10, and 11 we obtain 


% = (AyCr^Ai)-! + A[Q'^) (13) 


Equation 13 yields 

P = (AT Q-1 A^ )-i + (A[ Q-i A, )-l A{ Q-* A^Pj A^ Q"* Aj (AJ Q~^ A^ (14) 

Assume that the data noise is uneorrelated and that each data point has the same 
variance. Then 


Q = K 


( 15 ) 


where I is the identity matrix and is the common variance of the data. On the 
assumption that the values of the unadjusted gravity anomalies are independently 
distributed, the covariance matrix P 2 of gj can be written as 


o 


2 

1 


P 


2 


^^2 

0 


0 



(16) 


where n^ is the number of unadjusted gravity anomalies and a? is the second 
moment about zero of the i^*' unadjusted gravity anomaly. Also define a matrix 
K as 


K = (A[Q-iAj)-» aTQ-IAj (17) 

If is the number of adjusted parameters, then K is of dimension by Hj. 
With these assumptions Equation 14 yields the following expression for the vari- 
ance of the i*^ adjusted gravity anomaly 


«2 

* L ®ij »)>' 
j = 0 

where |3? ^ is the i*l^ diagonal element of the matrix (A| Aj )~^ and 

= K(I,j),j> 1 


(18) 


(19) 
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( 20 ) 


Define the error sensitivity matrix as 

S = 2. ••••. "pJ = 0. 1.2, 

And finally define the alias matrix as 

L = Sa (21) 

where 



The alias matrix reveals much of the probability structure of the estimation pro- 
cedure. From Equations 18, 20, 21 and 22 it can be seen that the standard devi- 
ation of the adjusted gravity anomaly is the root sum square of the terms in 
the i^^ row of the alias matrix. The elements in the first column of the alias 
matrix represent the R.S. S. contribution to the standard deviation of each esti- 
mated parameter due to the data noise. The elements in the jil* column, j ^ 2 
represent the R.S. S. contribution to the standard deviation of each estimated 
parameter due to the j - 1^* unadjusted parameter. These terms are called ^ 
aliasing contributions to the uncertainty in the adjusted parameters due to the 
uncertainty of the j - l*i unadjusted parameter. Notice that the aliasing con- 
tributions due to the j‘*^ unadjusted parameter are proportional to the standard 
deviai^pn of the j^*^ parameter. 

DEFINITION 

In a given estimation process the i^^ estimated parameter is 
said to be orthogonal with respect to the j*^^ unestimated 
parameter if the aliasing contribution to the i^l^ estimated 
parameter due to the uncertainty of the jll' unestimated param- 
eter is zero. It will be seen that the relation is symetric in 
the sense that if one parameter is orthogonal with respect to 
another and if they are interchanged within the adjusted and 
unadjusted modes, they are again orthogonal. 

To see the implications of the orthogonality relationship we need a more reveal- 
ing representation of the aliasing terms. Notice that the first term on the right 
side of Equation 14 is the covariance matrix of the estimation process under the 
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assumption that the unadjusted parameters are perfectly kno^. This covariance 
matrix gives the uncertainty of the estimates due only to the data noise. Define 
the so-called ’'noise only” covariance matrix as 


P = (A[ Q-i Aj)-i (23) 

Next observe that the elements in the i^*^ row ariH column of respectively Aj 
and A 2 are the partial derivatives of the i**' data point with respect to the ad- 
justed parameter and the partial derivative of the i^*^ data point with respect to 
the j**' unadjusted parameter. The aliasing contribution to the i^^ adjusted param- 
eter due to the j unadjusted parameter can be written as 


L(IJ + 1) 


k = 1 C = I 


3y(g) 

dg,(k) 


Q-'(C, 0 


3y(k) 

9l2(i) 


(24) 


vhere m is the number of data points. If the estimates of the adjusted param- 
eters are relatively uncorrelated in the noise only covairiance matrix, Equation 
24 can be approximated by 


L(IJ+ 1) 


m 


p(i.i) ^ 
e = 1 


3y(g) 

3g,(I) 


Q-'(C, C) 


3y(C) 

31-2 0 ) 


(25) 


From Equation 25 we see that the aliasing contribution is symetric in its argu- 
ments and this implies the symetry of the orthogonality relationship. 

A sufficient condition for the left side of Equation 25 to approximate zero is for 
the observability patterns of g, (I) and l’ 2 (j) in the data to be virtually nonover- 
lapping. Figures 1 and 2 demonstrate that the perturbation of a gradiometer 
observation due to a 1 mgal perturbation of a gravity anomaly rapidly attenuates 
with increasing spherical radius. Hence if the grids on which ^j(I) and g 2 (j) are 
defined are sufficiently separated, the orthogonality relationship would be effec- 
tively satisfied and the estimate of gj(I) would experience no aliasii^ from the 
uncertainty of g 2 (j). Conversely, if the grid on which gift) was defined were in 
close proximity to grids whose gravity anomalies were unadjusted, one would 
expect serious aliasing of the resultant estimate. 

It should be clear then, that if the gravity anomalies are estimated in a block the 
outer layers of the block contain gravity anomalies whose estimates will be badly 
aliased by the adjacent unadjusted parameters. It will be necessary to discard 
these estimates. But the gravity anomalies in a sufficiently small inner core of 
the block may be adequately separated from the unadjusted parameters to be ef- 
fectively orthogonal with respect to them. The estimates of these terms 
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GRADIOMETER OUTPUT PERTURBATION (IN ETVOS UNITS) X 10 



Figure 2. Gradiometer Output Perturbation Due to Imgal 
Perturbation of 2” Lunar Gravity Anomaly 
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presumably will be of sufficient accuracy that they can be accepted. In effect, 
for every block of gravity anomalies that we intend to estimate it will be neces- 
sary to construct of ’’buffer zone" several layers deep of gravity anomalies viiich 
surround the block. The new and larger block of gravity anomalies must be 
simultaneously estimated and then the estimates of gravity anomalies in the buf- 
fer zone must be rejected due to aliasing. To design an intelligent data process- 
ii^ procedure for a lunar gradiometer mission it is necessary to know the rela- 
tionship between the depth of the buffer zone and the accuracy of the estimation 
procedure. This relationship will vary with grid size, data block size, and satel- 
lite altitude. The most efficient way to study the relationship is to utilize covar- 
iance an^ysis techniques to generate alias matrices for several situations and 
to attempt generalizations from the results. This is accomplished in the next 
section. 
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RESULTS 


Assuming a satellite altitude of approximately 100 km and that a gradiometer 
reading is obtained every 15 seconds, a dense globally distributed set approxi- 
mately 86,000 gradiometer measurements are obtained mthin 15 days. In the 
low latitudes the sublunar points of the satellite at the measurement times form 
a grid on the moon of 75 degrees latitude by one degree longitude. All our simu- 
lations rely on this assumption concerning the data distribution. In doing so we 
eliminate the dynamical aspect of the problem and thus the aliasing effect of or- 
bit determination and attitude determination errors are not included in our 
restdts. 


As mentioned m the introduction the reason for eiujjlojnng gravity anom«uie8 to 
describe lunarpotential fine structure is the possibility of successfully estimat- 
ing local blocks of gravity anomalies in local blocks of gradiometer data. To 
investigate this possibility we focus attention on the estimation accuracy of a 
single gravity anomaly as a function of the number of adjacent gravity anomalies 
which are simultaneously adjusted and the quantity of data which is used in the 
estimation. We define the estimation block to be a block of gravity anomalies 
centered at the gravity anomaly in question and square in the sense that the block 
subtends the same number of degrees loi^tude as latitude. We assume that all 
the gravity anomalies in the estimation block are simultaneously adjusted in a 
given estimation procedure. As we have defined it, estimation blocks can only 
be incremented by integral numbers of twice the gravity anomaly block size. 

For 2® gravity anomalies, for Instance, possible estimation block sizes are 2°, 
6% 9% 10% 14% etc. 


The data block is defined as a gradiometer data set which was obtained when the 
sublunar points of the satellite were within a square region of the lunar surface 
centered at the gravity anomaly in question. For a given estimation block size 
the accuracy of a gravity anomaly estimate is not necessarily a monotone func- 
tion of the data block size. This may at first be puzzlii^ since in general one 
expects more data to yield better estimates. An explanation can be obtained by 
noticing that for a given estimation block the covariance matrix of the set of es- 
timated gravity anomalies is given by Equation 14 as the sum of a matrix which 
is dependent only on data uncertainty and a matrix which represents the aliasing 
effects from the unadjusted parameters. With increasii^ data blocks the ele- 
ments of the first matrix must decrease but in general tiie elements of the matrix 
which conveys the aliasing effects will increase. This effect can be shown gra- 
phically by means of so-called aliasii^ maps. To obtain the aliasing maps of 
Figures 3 and 4 we simulate in a covariance mode a situation in which we de- 
scribe lunarpotential fine structure by means of 2® by 2® gravity anomalies and 
we assume gradiometer data with an accuracy of .IE. In Figure 3 we mapped 
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Figure 4. Alias Map for 2° Gravity Anomalies 






































the R. S. S. contribution in mgals to the uncertainty in the estimate of the gravity 
anomaly defined on the blank grid element in the middle of the graph when the 
adjacent gravity anomalies are assumed to be in an unadjusted mode and to have 
an a priori uncertainty of 75 mgals. The data block as represented by the shaded 
area of the graph is 5® on a side. The satellite altitude was 50 km. Notice that 
die aliasing contributions decrease with the distance between an unadjusted param- 
eter and the estimated parameter thus demonstrating the inherent orthogonality 
properties of gravity anomalies in gradiometer data. In Figure 4 the data block 
size has been Increased to 20° and the aliasing effect of adjacent imadjusted grav- 
ity anomalies is seen to be greater. 

The choice of an intelligent data reduction strategy Is obviously dependent on a 
knowledge of the relationship between estimation accuracy and the choice of data 
block size and estimation block size. A computer program was written which 
computes the right side of Equation 14 for any given data reduction technique. 

The prc^ram was used to investigate data reduction strategies. First an alti- 
tude for the satellite was chosen and a gradiometer data block centered on a 
chosen gravity anomaly was assumed. The data accuracy was assumed as . IE. 

A normal matrix for 400 gravity anomalies was formed. The chosen gravity 
anomaly and gravity anomalies in successive layers surrounding the chosen 
gravity anomaly were assumed to be in an adjusted mode and the rest were placed 
in an unadjusted mode with uncertainties about zero of 75mgal. The right side 
of Equation 14 was computed and the resultant standard deviation in the estimate 
of the chosen gravity anomaly was recorded. The estimation block size was in- 
cremented by removing another layer of gravity anomalies from the unadjusted 
to the adjusted mode and repeatirg the process, Figxures 5 and 6 display the re- 
lationship between estimates of 1° and 2° gravity anomalies and the choices of 
data block size and estimation block size when the satellite is at an altitude of 
50 km. Figure 7 displays the information when the satellite is at 100 km. 

The graphs show that for 2° gravity anomalies an intelligent compromise be- 
tween computational load and accuracy is achieved when a data block size of 5° 
is employed and the buffer zone of rejected gravity anomaly estimates is at least 
two layers deep. The graphs show that when this strategy is followed and if the 
satellite is in a 50 km orbit, 2° gravity anomalies can be estimated with an ac- 
curacy of less than 1 mgal. For 1° gravity anomalies the best strategy utilizes 
a 10° data block size and a buffer zone of rejected gravity anomaly estimates 
five layers deep. With this strategy 1° gravity anomalies can be estimated to 
an accuracy of about 12 mgal. At an altitude of 100km, 1° gravity anomalies 
cannot be effectively estimated. Figure 7 shows that 2° gravity anomalies can 
be estimated from a satellite height of 100 km to an accuracy of 12 mgal provided 
a 20° data block size is utilized and provided that a buffer zone of rejected grav- 
ity anomalies five layers deep is utilized. 
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ASSUMPTIONS 

• A PRIORI ESTIMATES OF FREE 
AIR GRAVITY ANOMALIES 
±75 MILLIGALS 

• SATELLITE HT. = 50 KM 

• GRADIOMETER ACCURACY 0.1 E.U. 


Figure 5, Accuracy of 1* by 1® Mean Free Air Gravity Anomaly Estimates 
vs. Data Block Size for Various Estimation Block Sizes 
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Figure 6. Accuracy of 2® by 2^ Mean Free Air Gravity Anomaly Estimates 
vs. Data Block Size for Various Estimation Block Sizes 
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Figure 7. Accuracy of 2“ by 2° Mean Free Air Gravily Anomaly Estimates vs. Data Block 
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These results are dependent on the chosen value for the accuracy of the rotating 
gradiometer. Since this choice is somewhat arbitrary it is useful to parametric- 
ally represent gravity anomaly estimation accuracy as a function of gradiometer 
accuracy. Equation 18 implies that for a given choice of estimation block size 
and data block size, the standard deviation of a gravity anomaly estimate is the 
root sum square of a term which represents the aliasing effects of unadjusted 
gravity anomalies and a term ^ich is proportional to the standard deviation of 
the rotating gradiometer observation. Assuming a lO** data block size and that 
a given gravity anomaly is simultaneously estimated with all gravity anomalies 
in the five adjacent layers, the estimation accuracy for 1° gravity anomalies at 
an altitude of 50 km is given as a function of data accuracy as 

o(r, 50km) = (125440^ +0.342)^ (26a) 

At the same altitude and for 5° data block size and for a buffer zone of rejected 
gravity anomaly estimates two deep, accuracies of estimates of 2® gravity anoma- 
lies are given as a function of data accuracy as 

0(2°, 50km) = (90^ + 0.36)’^ (26b) 

At a height of 100 km and assuming a data block size of 20° and a buffer zone of 
rejected gravity anomalies five layers deep, the accuracies with which 2° gravity 
anomalies can be estimated are given as a function of data accuracy as 

0(2°, 100 km) = (16000 o2 + \p (26c) 

Standard deviations on the left side of the above equations are in mgal and the 
data standard deviation o on the right side is in Etvos units. 
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CONCLUDING REMARKS 


If a rotating gradiometer were placed on a satellite in a lunar polar orbit with an 
altitude in the 50 km to 100 km range, a dense and globally distributed set of 
gradiometer observations would be obtained in about 15 days. From this data set 
it should be possible to determine a global set of lunar gravity anomalies. But 
in order to keep the computational load within reasonable bounds it will be neces- 
sary to estimate local blocks of gravity anomalies in local blocks of gradiometer 
data. This report shows that if the satellite is flown at a 50 km altitude and if 
the proper estimation strategy is employed, an instrument with a . IE. accuracy 
permits estimates of l“ and 2° gravity anomalies which are accurate to within 
12 mgal and 1 mgal. When the satellite is flown at 100 km, 2° gravity anomalies 
can be estimated with an accuracy of 12 mgal. Equation 26 of this report per- 
mit a scaling of these results to a data accuracy different from .IE. 
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